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•  Explicit  algebraic  algorithm  which  can  model  any  charge-rest-discharge  protocol. 
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The  detailed  isothermal  electrochemical  model  for  a  lithium  ion  cell  has  ten  coupled  partial  differential 
equations  to  describe  the  cell  behavior.  In  an  earlier  publication  [Journal  of  Power  Sources,  222,  426 
(2013)],  a  reduced  order  model  (ROM)  was  developed  by  reducing  the  detailed  model  to  a  set  of  five 
linear  ordinary  differential  equations  and  nonlinear  algebraic  expressions,  using  uniform  reaction  rate, 
volume  averaging  and  profile  based  approximations.  An  arbitrary  current  profile,  involving  charge,  rest 
and  discharge,  is  broken  down  into  constant  current  and  linearly  varying  current  periods.  The  linearly 
varying  current  period  results  are  generic,  since  it  includes  the  constant  current  period  results  as  well. 
Hence,  the  linear  ordinary  differential  equations  in  ROM  are  solved  for  a  linearly  varying  current  period 
and  an  explicit  algebraic  algorithm  is  developed  for  lithium  ion  cell  voltage  prediction.  While  the  existing 
battery  management  system  (BMS)  algorithms  are  equivalent  circuit  based  and  ordinary  differential 
equations,  the  proposed  algorithm  is  an  explicit  algebraic  algorithm.  These  results  are  useful  to  develop  a 
BMS  algorithm  for  on-board  applications  in  electric  or  hybrid  vehicles,  smart  phones  etc.  This  algorithm 
is  simple  enough  for  a  spread-sheet  implementation  and  is  useful  for  rapid  analysis  of  laboratory  data. 

©  2013  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

An  isothermal  model  for  lithium  ion  cell  has  ten  coupled  partial 
differential  equations  (PDEs)  to  describe  the  cell  behavior  [1—4], 
namely:  two  solid  diffusion  equations  and  two  solid  current 
equations  in  the  electrode  regions;  three  electrolyte  diffusion 
equations  and  three  total  current  balances  in  the  negative  elec¬ 
trode,  separator  and  positive  electrode  regions,  respectively. 
Recently,  a  reduced  order  model  (ROM)  has  been  developed  [4], 
wherein,  using  uniform  reaction  rate,  volume  averaging  and  profile 
based  approximations,  the  aforementioned  ten  coupled  PDEs  are 
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reduced  to  a  set  of  five  linear  ordinary  differential  equations  (ODEs) 
and  nonlinear  algebraic  expressions.  The  nonlinearity  of  the  system 
once  relegated  to  algebraic  expressions  is  not  computationally 
cumbersome.  This  ROM  [4]  shows  negligible  errors  in  cell  voltage 
prediction  at  nominal  currents  and  less  than  3%  error  at  high  cur¬ 
rents  for  a  commercial  cell. 

The  current  work  derives  analytical  solutions  for  the  five  linear 
ODEs  in  the  ROM  [4],  thereby  developing  an  algebraic  algorithm  for 
cell  voltage  prediction  in  the  non-phase  change  lithium  ion  cells. 
Any  arbitrary  current  profile,  involving  charge,  rest  and  discharge, 
can  be  expressed  as  a  series  of  constant  current  and  linearly  varying 
current  periods.  The  latter  is  more  generic,  since  it  includes  the 
constant  current  period  as  well.  Hence,  analytical  solutions  are 
derived  for  a  linearly  varying  current  period  and  the  results  are 
developed  into  an  explicit  algebraic  algorithm  for  cell  voltage 
prediction.  This  algorithm  is  simple  enough  for  a  spread-sheet 
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implementation  and  can  be  used  in  battery  testing  equipment  for 
rapid  analysis  of  laboratory  data.  It  can  also  be  easily  developed  into 
a  battery  management  system  algorithm  for  on-board  applications 
in  electric  or  hybrid  vehicles  and  smart  phones. 

2.  Derivation  of  the  algebraic  ROM 


C,n(t)  =  Cln(ti)  -  ^  j/,(t  -  tj)  +  -  f*)2]  '  (10) 

Substituting  t  =  tf  and  simplifying  gives  the  average  concen¬ 
tration  of  solid  phase  Li  in  the  negative  electrode  at  the  final  time 


2.1.  Solid  phase  Li  concentration  and  state  of  charge  in  the  .  11/ 

electrodes  Cln(tf)  =  Cln(tj)  -  2pN^  (tj  +  /f)  (tf  -  tj).  (11) 


The  evolution  of  average  solid  phase  Li  concentration  in  the 
negative  electrode  (Eq.  (127)  in  Ref.  [4])  is 


d  cln  3(?n) 

dt  rn  ' 


(1) 


The  average  reaction  rate  in  the  above  equation  is  related  to  the 
external  current  and  the  electrode  design  parameters  as  (Eq.  (14)  in 
Ref.  [4]) 


<jn)(t)  = 


m 

anFln  ’ 


(2) 


Using  the  condition,  Ii  =  I{  =  I  the  above  expression  reduces  to 
the  solution  for  a  constant  current  period. 

According  to  the  electro-neutrality  condition,  for  every  neutral 
Li  ionized  in  an  electrode,  a  Li  +  ion  is  neutralized  in  the  opposite 
electrode.  Thereby,  the  total  number  of  neutral  Li  in  the  solid  phase 
of  both  the  electrodes  remains  constant  (Eq.  (128)  in  Ref.  [4]) 

/neinCin(f)  +  ipcipCip(t)  =  /nt’uiCino  +  (pcipCipo  =  Constant. 

(12) 

In  terms  of  scaled  quantities,  the  above  equation  becomes 


where,  an  is  the  specific  surface  area  of  the  active  material  in 
negative  electrode  given  by, 


Using  Eqs.  (2)  and  (3)  in  (1)  gives 


dt  IneinF 


(4) 


Defining  the  scaled  average  solid  Li  concentration  in  negative 
electrode  as 


and  the  maximum  specific  moles  of  solid  phase  Li  in  the  negative 
electrode  as 


NlnCm(t)  +  NlpClp(t)  =  NlnSoCn0  +  NipSoCpo,  (13) 

where  Nlp  =  /pcipcsp  max-  (14) 

The  initial  condition  is  assumed  to  be  a  completely  equilibrated 
state,  that  is,  no  concentration  gradients  or  equivalently,  the 
average  concentration  of  Li  within  the  active  material  spheres  is 
equal  to  the  surface  concentration.  Hence,  the  initial  scaled  con¬ 
centrations  are  identical  to  initial  electrode  states  of  charge.  Using 
the  Eq.  (13),  the  average  Li  concentration  in  the  positive  electrode  is 
computed  as 


Cip(f)  =  SoCpo  +  ^  [socn0  -  Cln(f)] .  (15) 

The  average  concentration  gradient  within  the  active  material 
spheres  in  the  negative  electrode  is  scaled  using  maximum  con¬ 
centration  and  particle  radius  as 


Nln  =  /neinCsnmax,  (6) 

Eq.  (4)  is  rewritten  in  terms  of  scaled  quantities  as 


dCln  J(t) 
dt  FNln ' 


(7) 


Time  is  left  unsealed  since  there  are  many  different  time  scales 
in  this  problem.  Consider  a  linear  change  of  external  current  from  L 
to  If,  in  the  time  interval  tj  to  tf.  The  current  at  any  time  in  this  in¬ 
terval  is 


l(t)  =  4  +  7 - f(t-ti)- 


(8) 


The  evolution  equation  for  the  radial  concentration  gradient  in 
negative  electrode  is  given  by  (Eq.  (153)  in  Ref.  [4]) 


^  +  ^Clm  =  -^(,n>.  (17) 

With  the  solid  diffusion  time  scale  in  negative  electrode  defined 


Using  Eqs.  (2),  (8),  (16)  and  (18), 


Using  the  above  expression  in  Eq.  (7)  and  integrating  from  tj  to  t 
gives 

Cin(t)-Cln(tj)  =  -FA[|n  |fj  J dt+^^y  J (t-tj)dtj  .  (9) 

On  integrating  and  rearranging  the  above  expression,  the  scaled 
average  concentration  of  solid  phase  Li  in  the  negative  electrode  at 
any  time  is 


dCirn  30  ^ 
dt  +TlnClrn 


15  1 
T  FN,n 


M) 


(19) 


Eq.  (19)  is  a  first  order  differential  equation  which  can  be  solved 
by  more  than  one  method.  M.  R.  Spiegel  [5]  has  details  of  such 
solution  methods.  In  the  present  work,  the  equation  is  solved  by 
obtaining  the  general  solution  which  is  the  sum  of  a  complemen¬ 
tary  function  and  a  particular  integral.  The  complementary  func¬ 
tion  obtained  by  solving  the  homogeneous  ODE,  is 


Cirn  =  anexp  -30  — 
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(20) 


!p  +  ?p(t-ti), 


To  obtain  the  particular  integral,  a  trial  function  linear  in  time  is 
substituted  in  Eq.  (19) 

Cl'm  =  0n  +  Yn(t-ti)-  (21) 

Equating  the  coefficients  of  the  constant  and  the  linear  terms 
in  time,  in  Eq.  (21),  the  expressions  for  /3n  and  yn  are  obtained 


l(lf-li)tip  i 

where  y0  =  — — - - - - 


-  TP 
_  1  A^ip  1 


FNip 


iH 


jTptip, 

and  «p  =  C|rp(tj)  — /Ip. 


(34) 

(35) 

(36) 


,  O  1  fitjn  1 

3nd  Pn  = 


The  complete  solution  is  the  sum  of  complementary  function 
and  particular  integral 

Cim(t)  =  an  exp  [  -  30 (f  ^1  +  /?n  +  y n(t  -  t{).  (24) 

L  Tln  J 

At  the  initial  time,  the  above  expression  reduces  to 
Clm(ti)  =  <*n  +  Pn,  (25) 

giving  an  =  Clm(tj)  -  /Jn.  (26) 

Thus,  Eqs.  (22),  (23)  and  (26)  determine  the  three  coefficients  in 
the  complete  solution  of  Eq.  (24).  For  the  initial  condition  of  an 
equilibrated  state,  the  initial  average  radial  concentration  gradient 
can  be  set  to  zero.  Knowing  the  average  concentration  and  the 
average  radial  concentration  gradient,  the  negative  electrode  SoC  is 
obtained  by  scaling  Eq.  (149)  in  Ref.  [4]  as, 


Using  the  expressions  for  the  average  concentration  and  the 
average  radial  concentration  gradient,  the  positive  electrode  SoC  is 
obtained  by  scaling  Eq.  (149)  in  Ref.  [4]  as 

SoW  ’  SS  *  +BC'-(t)  +  T®T|r  (37> 


2.2.  Electrolyte  phase  Li  concentration  in  the  electrodes  and  the 
separator 


The  initial  specific  moles  of  Li,  in  the  electrolyte  phase  of  the 
electrode  and  separator  regions,  are  defined  as 


and  the  electrolyte  diffusion  time  constants  x2k  =  (39) 

D2k 

where  k  =  n,s,  p.  It  is  to  be  noted  that  in  Eq.  (39),  D2k  is  the  effective 
diffusivity  that  includes  the  effects  of  porosity  and  the  Bruggeman 
factor.  Details  of  the  scaling  are  given  in  the  ‘Notation’  section.  The 
total  specific  moles  of  Li  in  the  electrolyte  phase  of  the  cell  is  given  by 


N2  =  N2n+N2s+N2  p. 


(40) 


SoCn(t)  m 


Csn(t) 


=  Cln(t)+_Clnr(t)-w 


1  /(t)Tln 

105  FNln  ' 


(27) 


Similarly,  the  average  radial  gradient  in  the  positive  electrode  is 
obtained  by  solving  the  following  equations  (Eqs.  (153)  and  (15)  in 
Ref.  [4]) 


Electrolyte  concentration  field  is  characterized  by  the  electro¬ 
lyte  flux  at  the  two  electrode-separator  interfaces:  q2in(t)  and 
<?2ip(t).  The  coupled  ODEs  for  these  two  internal  variables  in  ROM 
are  (Eqs.  (83)  and  (84)  in  Refs.  [4],  along  with  the  average  reaction 
rates  Eqs.  (2)  and  (5)  herein) 


dc1rp  30P1p  _  45/  X 

'  dt "+  r*  Clrp  “  2r^vp/’ 

(28) 

(41) 

where(Jp)(t)  - 

(29) 

and  a2 1cl^  +  =  <?2iP  -  (1  -  4>^r  ’ 

(42) 

where  the  coefficients  (having  units  of  time)  in  these  ODEs  are 

3ei„ 

written  in  terms  of  the  scaled  quantities  as 

and  ap  =  — p . 

(30) 

rp 

„  (1  N2n  t2s  r2n\  Nm  ,  1  N2n  *2n 

11  “  2s  ~6~  ^7  N^+2A^T2s  +  X’ 

(43) 

With  the  maximum  solid  phase  Li  capacity  from  Eq.  ( 14)  and  the 

solid  diffusion  time  constant 

rp 

Tlp  =  iV 

(31) 

„  (1N2n  t2s  r2p\  N2n  LN2n 

12  \2N2s  2s  +  3  3  )  N2  +2N2s  2s’ 

(44) 

Eq.  (28)  gets  scaled  as 

n  f1N2nT  t2s  r2n\  N2p 

21  =  \2 Ws  2s+~e  W* 

(45) 

(32) 

and  a22  =  (' N2n  ~2s  ■  T2s  Nlp  t2p  . 

22  \2  N2s  2s  3  3  )  N2  3 

(46) 

Its  complete  solution  is  obtained  as 

Rearranging  Eq.  (41 )  gives 
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d(?2ip 

dt 


;  n^q2in  o' 


an  dq2in  |  1  (1  -jj 

at2  dt  a12  F 


Using  this  equation  in  Eq.  (42)  and  rearranging  gives 


Differentiating  the  above  equation  with  respect  to  tir 
using  the  current  profile  of  Eq.  (8)  gives 


(48) 


d(?2ip  a22  dq2in  detQd2q2in  /a^  \  (1  - 1+)  //f-fA 

dt  aX2  dt  a12  dt2  \ai*  /  F 

(49) 

Using  the  above  equation  coupled  with  Eq.  (8)  in  Eq.  (41)  fol¬ 
lowed  by  simplification  gives 


y(t-tj-a12-a22)  • 


(50) 


This  equation  is  a  second  order  ODE  with  a  time-dependent 
forcing  function.  Its  complementary  function  is  obtained  from  the 
homogeneous  part  of  the  ODE  as  follows:  Solving  the  quadratic 
equation 


detaA2  -  (an  -  a22)A  -1=0, 


(51) 


gives  the  roots  Ax  and  A2.  Then,  the  complementary  function  is 


determined  from  the  initial  value  and  the  initial  rate  of  change  of 
q2in,  as  shown  below.  Differentiating  Eq.  (57)  gives 

=  AAX  exp(Ajt)  +  BA2  exp(A2t)  +  Q_.  (58) 

Using  Eqs.  (57)  and  (58)  in  Eq.  (48)  and  simplifying  gives 

<72iP  =  -  1  (a22  +  detaA,  )A  exp(A,  t)  -  1  (a22 
H  ai2  ai2 

+  detaA2)B  exp(A2t)  -  —  (P  +  Qt)  -  Q+  (— 
a\2  a  u  \ai2 

+  i)  ^1(0. 

(59) 

The  initial  values  of  q2jn  and  q2,p  ate  set  to  zero  at  the  beginning 
of  the  simulation  based  on  the  assumption  of  equilibrated  state.  At 
successive  time  intervals  they  are  obtained  by  their  continuity  in 
time  i.e.  their  values  at  the  end  of  the  previous  time  interval  are  the 
initial  values  for  the  successive  time  interval.  Once  the  initial  values 
of  q2in  and  q2,p  are  known,  the  coefficients  A  and  B  can  be  obtained 
by  solving  Eqs.  (57)  and  (59). 

3.  Algebraic  ROM  algorithm 

In  this  section  an  explicit  algebraic  reduced  order  algorithm  is 
constructed  using  the  theoretical  results  collated  in  Section  7  of  [4] 
and  those  derived  in  the  earlier  section  of  this  work.  The  initial  state 
is  assumed  to  be  equilibrated  i.e.  the  cell  is  rested  until  the  con¬ 
centration  gradients  in  the  electrolyte  and  solid  phases  have 
relaxed.  The  steps  in  the  algorithm  are: 


q2in(t)  =  Aexp(Ait) +  Bexp(A2t).  (52) 

The  particular  integral  is  obtained  as  follows:  Substituting  the 
linear  trial  function 

qSn(0  =P  +  Qt,  (53) 

in  the  ODE  Eq.  (50)  and  simplifying  gives 

(an-a22)Q+(P  +  Qt)  =  (k_4)  L  +  (/f  Jl)  (t_t._ai2 

f  L  (tf-ti) 

-  a2 2)J . 

(54) 

Equating  the  coefficients  of  the  linear  and  constant  terms  in 
time  gives 

a_g-4)('f-',)  (55) 

F  ("-'■) 

and,  P  =  (1  ~4)fi  -  (ti  +  au  +  a12)Q.  (56) 

The  complete  solution  to  the  ODE  Eq.  (50)  is  the  sum  of  the 
complementary  function  Eq.  (52)  and  the  particular  integral  Eq. 

(53) 

<72in(t)  =  A  exp(Ax  t)  +  B  exp(A2t)  +  P  +  Qt.  (57) 

With  P  and  Q.  determined  from  Eqs.  (55)  and  (56),  A  and  B  are 


1.  Computation  of  average  concentration  based  SoC  in  negative 
electrode:  Eq.  (10). 

2.  Computation  of  average  radial  gradient  in  negative  elec¬ 
trode:  Eq.  (24). 

3.  Computation  of  SoC  in  negative  electrode:  Eq.  (27). 

4.  Computation  of  average  concentration  based  SoC  in  positive 
electrode  using  solid  Li  balance:  Eq.  (15). 

5.  Computation  of  average  radial  gradient  in  positive  electrode: 
Eq.  (33). 

6.  Computation  of  SoC  in  positive  electrode:  Eq.  (37). 

7.  Computation  of  interfacial  electrolyte  fluxes:  Eqs.  (57)  and 
(59). 

8.  Computation  of  interfacial,  average,  collector-end  and  mid¬ 
separator  electrolyte  concentrations  using  interfacial 
fluxes:  Eqs.  (77),  (78),  (57),  (74),  (65),  (60),  (67),  (76)  of  [41. 

9.  Computation  of  interfacial  and  collector-end  electrolyte  po¬ 
tentials:  Eqs.  (105),  (106),  (114)  and  (117)  of  [4], 

10.  Computation  of  collector-end  solid  potentials  and  hence  the 
cell  voltage:  Eqs.  (157)  and  (158)  of  [4], 

This  algorithm  is  used  to  compare  the  algebraic  ROM  prediction 
with  the  Matlab-Simulink  result  of  the  original  ROM  of  [4]  for  a 
constant  current  discharge  case  in  Fig.  1.  This  figure  shows  that  the 
order  reduction  to  algebraic  level  is  indeed  exact  at  both  nominal 
and  high  current  discharge  rates. 

4.  Conclusions 

Real-time  voltage  prediction  of  lithium  ion  cells  using  the 
electrochemical  model  is  difficult  owing  to  solving  coupled  partial 
differential  equations.  In  an  earlier  publication  [4]  a  reduced  order 
model  (ROM)  was  developed  by  reducing  the  partial  differential 
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0  500  1000  1500  2000  2500  3000  3500 

Time  (s) 

Fig.  1.  Comparison  of  Excel  implementation  of  algebraic  ROM  with  the  Simulink 
implementation  of  ODEs  based  ROM  for  constant  current  discharge  at  various 


equations  to  a  set  of  five  linear  ordinary  differential  equations  and 
nonlinear  algebraic  expressions.  This  is  achieved  [4]  by  applying 
order  reduction  methodologies  such  as  uniform  reaction  rate 
approximation,  volume  averaging  and  profile  based  approxima¬ 
tions.  In  the  present  work,  this  model  reduction  approach  is  applied 
for  any  arbitrary  current  profile.  An  arbitrary  current  profile, 
involving  charge,  rest  and  discharge,  is  expressed  as  a  series  of 
constant  current  and  linearly  varying  current  periods.  The  latter  is 
more  generic,  since  it  includes  the  constant  current  period  as  well. 
The  linear  ordinary  differential  equations  in  ROM  are  solved  using 
the  analytical  form  of  the  linearly  varying  current  and  an  explicit 
algebraic  algorithm  is  developed  for  lithium  ion  cell  voltage  pre¬ 
diction.  Owing  to  the  simplicity  of  the  algorithm,  it  is  amenable  for 
implementation  on  spreadsheets  which  can  be  used  to  validate 
laboratory  test  data  and  can  be  extended  to  a  battery  management 
system  for  onboard  application  in  electric  vehicles  and 
smartphones. 
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Notation 


Constants  &  variables 

R:  Gas  constant,  8.314  J  mol-1  K-1 

F:  Faraday's  constant,  96487  C  mol-1. 

T:  Temperature,  298  K  in  this  work. 

I  (t):  External  current  density  at  time  t,  A  m-2. 

If.  External  current  densities  at  times  q  and  tf  in  a  linearly  varying  current  pulse, 

t+:  Electrolyte  transference  number,  dimensionless. 

brug:  Bruggeman  factor,  used  to  calculate  effective  porous  media  properties, 
dimensionless. 

l„,ljp:  Thickness  of  negative  electrode,  separator  and  positive  electrode,  m. 
L  =  ln  +  ls  +  lp  is  total  cell  thickness. 

r„,  rp:  Radii  of  active  material  spheres  in  the  negative  and  positive  electrodes,  m. 
eik,e2k,ejk:  Volume  fractions  of  active  material,  electrolyte  and  filler  in  electrode  or 
separator  domains,  dimension  lessen,  +  c2k  I  ffk  =  L  Separator  has  no  active 
material.  Hence  eis  =  0. 

ak  =  3ejk/rk:  Specific  surface  area  of  active  materials  in  n  and  p  electrodes,  m-1. 

jk:  Local  surface  reaction  rate  in  k  =  n,p  electrodes,  mol  m-2s-1. 

iko  =  kktekjiax  -  csk)a5c3fc25-'  Over-potential  independent  rate  pre-factor, 

kk:  Surface  reaction  rate  constant  in  k  =  n,p  electrodes,  (mol  m-2  s-1)(mol  m-3)-1'5. 
Ci,  c2:  Solid  phase  concentration  in  active  material  spheres  and  electrolyte  concen- 

csk,  Csk  max ■'  Surface  solid  phase  concentrations,  and  their  maximum  value,  mol  m-3. 
Ckk  =  cik/csk  max'  Scaled  average  solid  phase  concentrations  in  k  =  n,p  electrodes, 
dimensionless. 

C\rk  =  tkCirt/Cstmax-'  Scaled  average  radial  gradient  in  k  =  n,  p  electrodes, 
dimensionless. 

C2ikA2ik-  Interfacial  electrolyte  concentraions  and  fluxes,  mol  m3,  mol  m-2  s-1. 

D2:  Electrolyte  diffusivity  (material  property),  m2  s-1. 

D2 k  =  P>2  x  fj™8:  Effective  electrolyte  diffusivity  in  the  porous  k  =  n,s,p  domains. 
Dik:  Solid  diffusivity  of  k  =  n,p  electrode  material  spheres,  m2  s-1. 
k2:  Electrolyte  conductivity  (material  property),  S  m-1. 

K2k  =  *2  x  e2™S;  Effective  electrolyte  conductivity  in  porous  k  =  n,  s,  p  domains. 
Nlk:  Maximum  specific  moles  of  solid  Li  in  k  =  n,p  electrodes,  mol  m-2. 

N2k:  Initial  specific  moles  of  electrolyte  Li  in  k  =  n,  s,  p  regions,  mol  m-2. 

N2  =  JV2„  +  N2s  +  N2p:  Total  molar  capacity  of  electrolyte  Li  in  the  cell,  mol  m-2. 

<rk:  Electrical  conductivity  of  k  =  n,p  electrode  materials  (material  property),  S  m-1 
<tjk  =  ak  x  C]™8.'  Effective  electrical  conductivity  in  the  porous  k  =  n,  p  domains, 

4>i,  $2:  Solid  and  liquid  potential,  V. 

rjk:  Solid  diffusion  time  constants  in  k  =  n,p  electrodes. 

1 2k:  Electrolyte  diffusion  time  constants  in  k  =  n,s,p  domains. 

Uk:  Open  circuit  voltage  of  k  =  n,p  electrodes,  V. 

SoCk  =  Cjjj/c,,.  mox  k:  State  of  charge  of  k  =  n,p  electrodes,  dimensionless. 

Subscripts  and  operators: 

1,  2,  f:  Solid  active  material,  liquid  electrolyte  phases  and  solid  fillers  (including 
additives) 

k  =  n,s,  p:  Negative  electrode,  separator,  and  positive  electrode. 
in:  Separator  interface  of  the  negative  electrode. 
ip:  Separator  interface  of  the  positive  electrode. 

CF:  Complementary  function  of  an  ordinary  differential  equation. 

PI:  Particular  integral  of  an  ordinary  differential  equation. 

0:  Initial  condition. 

(•):  Volume  average,  defined  in  the  three  domains,  n,  s.  p. 

A  Volume  average  within  an  active  material  sphere  in  n  or  p  electrode. 


